Are Percolation Transitions always Sharpened by Making Networks Interdependent? 
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We study a model for coupled networks introduced recently by Buldyrev et al, Nature 464, 
1025 (2010), where each node has to be connected to others via two types of links to be viable. 
Removing a critical fraction of nodes leads to a percolation transition that has been claimed to be 
more abrupt than that for uncoupled networks. Indeed, it was found to be discontinuous in all cases 
studied. Using an efficient new algorithm we verify that the transition is discontinuous for coupled 
Erdos-Renyi networks, but find it to be continuous for fully interdependent diluted lattices. In 2 
and 3 dimension, the order parameter exponent j3 is larger than in ordinary percolation, showing 
that the transition is less sharp, i.e. further from discontinuity, than for isolated networks. Possible 
consequences for spatially embedded networks are discussed. 
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While the theoretical study of single networks has ex- 
ploded during the last years, relatively little work has 
been devoted to the study of interdependent networks. 
This is in stark contrast to the abundance of coupled net- 
works in nature and technology - one might e.g. think of 
people connected by telephone calls, by roads, by their 
work relationships, etc. For single networks it is well 
known that removing nodes can lead to cascades where 
other nodes become dysfunctional too and deleting 
a sufficient fraction of nodes leads to the disappearance 
of the giant connected cluster. If the network is already 
close to the transition point, deleting a single node can 
lead to an infinite cascade similar to the outbreak of a 
large epidemic in a population. 

Assume now that all nodes have to be connected via 
different types of links in order to remain functional. It 
was argued in Q that in such cases the cascades of fail- 
ure triggered by removing single nodes should be greatly 
enhanced, and that the transition between existence and 
non-existence of a giant cluster of functional nodes should 
become discontinuous. This claim was backed by a mean 
field theory that becomes exact for locally tree-like net- 
works (e.g. large sparse Erdos-Renyi (ER) networks), 
and by numerical simulations for various types of network 
topologies. In the present paper we show that this view 
is not entirely correct: For fully interdependent diluted 
rf-dimensional lattices, the transition is not only contin- 
uous, but it is less sudden than the ordinary percolation 
(OP) transition for isolated lattices and represents a new 
universality class. 

The problem is best illustrated by an actual case dis- 
cussed in H, which concerns an electric power blackout 
in Italy in September 2003 Q • According to Q (see also 
[3, [1|), the event was possibly triggered by the failure of a 
single node iq in the electricity network. Nodes in a power 
networks are in general also linked by a telecommunica- 
tion network (TN) and need to receive information about 
the status of the other nodes. In the present case, pre- 
sumably some nodes in the TN failed, because they were 



not supplied with power. This then led to the failure 
of more power stations because they did not receive the 
necessary information from ig, of more TN nodes because 
they were not supplied with electric power, etc. The en- 
suing cascade finally affected the entire power grid. 

The crucial point here is that each node has to be con- 
nected to two distinct networks that provide different ser- 
vices, in order to be viable. At the same time nodes act 
as bridges to bring supply to other nodes. If a node gets 
disconnected from one network, it no longer can func- 
tion and looses also its ability to serve as a connector 
in the other. The claim in to be scrutinized here, is 
that these cascades of failure are much more abrupt in 
interdependent networks than in isolated ones, leading to 
much sharper transitions. 

In a single network, the existence of an "infinite" clus- 
ter of nodes, making possible the outbreak of a large 
epidemic, is described by OP. Whether such a large out- 
break can happen depends on the average connectivity 
of the network, characterized by some parameter p. If p 
is below a critical value p c , no infinite epidemic can oc- 
cur, while it occurs with probability P > if p > p c . 
For p slightly above p c , both P and the relative size 
of the epidemic in a large but finite population scale 
~ (p — Pc) , where the order parameter exponent j3 de- 
pends on the topology of the network. For ER networks 
P = 1, while for randomly diluted d-dimensional lattices 
P depends on d, with /?(d = 2) = 5/36 « 0.1389 @ and 
j3{d = 3) = 0.4170(3) 0]. In all these cases /3 > 0, mean- 
ing that the transition is continuous. A discontinuous 
transition, as found in would correspond to fj = 0. 

Discontinuous percolation transitions have recently 
been claimed to exist in several other models H 0] , in- 
cluding explosive percolation [|| . The numerical evidence 
for discontinuity given in Q was supported in numerous 
papers. It became clear only recently that the transition 
is actually continuous, although with small (3 and with 
unusual finite size behavior [10( . In view of the difficulty 
to distinguish numerically between a truly discontinuous 
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transition and a continuous one with very small /3, we 
decided to perform more precise simulations. 

The algorithm used in [2| follows in detail the cascades 
triggered by removing nodes and, as a result, does not 
allow one to study large networks with high statistics. In 
our simulations, instead of removing nodes, we add nodes 
one by one. Using a modification of the fast Newman- 
Ziff algorithm 
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] , this gives a code which no longer fol- 
lows entire cascades, as they are broken up into short 
sub-cascades, and gluing them together would make the 
algorithm slow again. But it allowed us to obtain high 
statistics for reasonably large systems. 

The model is formally defined as follows: Start with a 
single set N of N nodes and with two networks A and B 
that are obtained by linking these nodes (notice that A 
and B need not be connected, and indeed some nodes in 
M may be not connected at all, in which case A and B 
make use only of subsets of W; also we do not demand 
that all links in A and B are different). Typically, we 
construct A and B by starting with a dense network and 
deleting randomly links from it, keeping links only with 
probability q < 1. In this way, ER networks are con- 
structed by starting with a complete graph and keeping 
only L = qN(N — l)/2 links. Alternatively, diluted reg- 
ular d-dimensional lattices are obtained by starting with 
a (hypcr-)cubic lattice with N = L d nodes and helical 
boundary conditions, and keeping only a fraction q of 
the dN links. 

On these coupled networks (each obtained by bond 
percolation with parameter q), we study a site per- 
colation problem by retaining only a fraction p of all 
nodes, calling the set of retained nodes Af p . We de- 
fine A6-clustcrs as subsets of nodes £ Af p that are con- 
nected both in A and in B. More precisely, assume that 
C = i-2, ■ ■ ■ 4} is a subset of nodes in J\f p . We call it 
a (connected AB-) cluster, if any two points i £ C and 
j £ C are connected by (at least) two paths: one path us- 
ing only links G A, and nodes only £ C, and another path 
using only links e B, also using nodes only £ C. Notice 
that we do not allow paths that involve nodes outside C, 
i.e. *4i3-clusters are 'self-sustaining'. The "order param- 
eter" S = rrimnx/N is then the relative size of the largest 
>4£>-clustcr, for given p and q. 

To find these maximal clusters, we start with an empty 
initial configuration with no nodes but with a list of all 
possible links in A and£>, and set m max = 0. Then we 
add nodes one by one. Each time a new node i is added, 

(a) We check whether it is linked to any of the existing 
nodes. If it is not linked to any other node either by A 
or by B links, we simply insert the next node. 

(b) Otherwise, we update the cluster structures in A 
and B separately by means of the Newman-Ziff algorithm, 
and denote the sets of nodes linked to i by Ca an d Cg. If 
one of them has size < m max , then m max cannot increase 
and we insert the next node. 

If not, we check whether the biggest „4Z?-cluster in 



Ca D can have a size > m max , by following a cascade 
similar to that in 0]. If the cascade stops at a cluster 
size > m max , then m max is increased. If it continues to 
a size < ?ri max , the cascade is stopped and m maK is left 
unchanged. In cither case, we then insert the next node. 

(c) This process continues until a preset value p niax is 
reached. Stopping at p < 1 is crucial for efficiency, as 
the algorithm slows down dramatically at large p. We 
typically follow the evolution up to p slightly above p c 
for all realizations, and follow it up to larger values of 
p for successively fewer runs. This reflects the fact that 
simulations are slow for p^> p c , but fluctuations are also 
smaller, so that fewer samples arc sufficient. 

For ER graphs the model can be simplified, since bond 
and site dilution both lead again to ER graphs. Hence we 
do not have to distinguish between them and can skip the 
site percolation part. The order parameter S = 77i max /./V 
is then, in the limit N — > oo, a unique function of the 
average degree (k). This function is easily found by ar- 
guments analogous to those for single networks. 

Consider an isolated ER network with average degree 
(k) = z in the regime where an infinite cluster exists, i.e. 
where an infection has a non-zero chance to lead to an 
infinite epidemic. Let Si be the probability that node i 
gets infected during this epidemic. The probability that 
i docs not get infected is then 



i - St = n a - s'j), 



(i) 



where the product runs over all neighbors of i. Here, 
S'j is the probability that j is infected, conditioned on 
it being picked as a node at the end of a link, and we 
used the fact that the graph is locally tree-like, so all S'j 
are independent. For ER graphs the degree distribution 
is Poisson, and S' and S obey the same statistics. Av- 
eraging Eq. |T|) over all nodes and topologies gives then 



12. 



= y^(i 



ki 



S) k 



-zS 



(2) 



where we dropped the index on S. Otherwise said, the 
probability S that any site is linked to the infinite cluster 
is 1 — exp(— zS). For two interdependent ER networks 
with average degrees za and zb, the chance to belong to 
the infinite „4£>-cluster is equal to the probability to be 
linked to it both via A and via B, giving 



S= (1 -e~ ZAS ){l -e- ZBU ) 



—zbS\ 



(3) 



Although this is much simpler than the theory presented 
in Q, it is exactly equivalent. It is generalized trivially 
to > 2 interdependent networks [lij], and to other types 
of interdependencies (H)]. If z_a = zb = z, one finds 
only the solution S = for z < z c = 2.455407 . . ., while 
a second stable solution S > exists for z > z c . Just 
above threshold, S c = 0.511699 . . .. 
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FIG. 1. (Color online) Plot for S = (m m ax)/iV against z, for 
two interdependent ER networks with degrees za — zb = z. 
For technical reasons, each curve does not correspond to a 
fixed value of N, but of No = 4JV/jz. The grey curve is the 
solution of Eq. Q. The intersection of the horizontal and 
vertical lines indicate the point (z c ,S c ). 
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FIG. 2. (Color online) Scatter plot for m max /iV against z, 
just after the largest jump in m max . Color corresponds to a 
fixed value of N. The lines indicate the analytic prediction 
for the point (z c ,S c ), according to Eq. ([3]). 

Results from our numerical simulations for ER graphs, 
using the algorithm outlined above, are shown in Figs. [T] 
and [21 Figure [1] shows S versus z for networks of differ- 
ent sizes. Each curve is based on 10 4 runs, except for 
the largest N. The data indeed approach the theoreti- 
cal curve (indicated in grey), as N — > oo. While Fig. Q] 
demonstrates that the theory gives the correct z c , it is 
much harder to argue that it gives also the correct S c . 
To see this, we notice that m max /N makes in each run 
exactly one big jump, from rss to rss S c - The values of z 
and S just after the jump are shown as scatter plots in 
Fig. [2] We see clouds of points that are indeed centered 
near z c and S c , and whose sizes decrease with N. 

For bond percolation on the square lattice, the OP 
threshold is at q c = 1/2 Q. We therefore look for ^IB- 
percolation in the parameter range 1/2 < q < 1. We 
assume the usual finite size scaling (FSS) ansatz Q 

(m max )=L D ff[(p-p c )L 1 /»], (4) 

where v is the correlation length exponent, Df = d — /3/i> 
is the fractal dimension of the incipient infinite cluster, 
and f{z) is a smooth (indeed analytic) function. Accord- 



FIG. 3. (Color online) Data collapse for (m ma x)/i I3/ against 
(p— p c )I/ 1//y , for 2-d lattices. Each set of curves corresponds to 
one value of q, while each curve within each set corresponds to 
a system size L. For this plot, v = 1.19 and Df = 1.85 were 
used, and the values of p c are 0.96025, 0.77556, and 0.6544 
for q = 0.6, 0.75, and 0.9. Due to universality of the scaling 
function f(z) we can also collapse the three set of curves, 
by multiplying m max /-£/ D/ and (p — p^L 1 '" by suitable in- 
dependent factors (data not shown). 
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FIG. 4. (Color online) Critical values p c versus q for two 
coupled 2-d lattices. Error bars are much smaller than the 
symbol sizes. Notice that the bond percolation threshold on 
the square lattice is q c — 1/2, thus the curve cannot extend 
below q < 0.5. The transition is in the same universality class 
even for q — q m ; n = 0.5757(4), where p c — 1 and the model 
simplifies, as no site percolation is involved. For q — ¥ 1, the 
model crosses over to OP. 

ing to this ansatz, we expect a data collapse if we plot 
(m max )/L z;i/ against (p — p c )L x ^. Three such data col- 
lapses are shown in Fig. [3l each for a different value of 
q. Each of the three "curves" in this figure are indeed 
several collapsed curves corresponding to different val- 
ues of L in the range 2 5 to 2 9 , obtained from more than 
10 6 realizations for the smallest lattice and w 10 4 for the 
largest. For all curves the same values of Df and v were 
used, while p c depends of course on q. The values of p c 
are plotted against q in Fig. 0J 

The fact that data collapse was obtained in Fig.[3]foi' g- 
independent values of the exponents indicates that these 
exponents are universal for g m j n < q < 1. But a closer 
inspection of Fig. [3] shows that the quality of the collapse 
deteriorates as q — > 1, due to the expected cross-over 
to OP (for q —¥ 1, A and B become identical, and the 
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FIG. 5. (Color online) Log-log plot of {mma.x)/L Df against 
L, for 2-d interdependent percolation at q — 0.60 and at fixed 
values of p. At the critical point (p c = 0.96025(20)) we expect 
a straight line. The value of the fractal dimension Df is chosen 
such that this line is horizontal. 
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FIG. 6. (Color online) Log- log plot of S = (m max )/L against 
p — pc, for 2-d interdependent percolation with q = 0.60. The 
slight upward curvature for large p — p c indicates the limit of 
the critical region, while the upward curvature for p — p c — > 
is due to finite size corrections. 



problem crosses over to OP). Thus we use data for q = 0.6 
for more detailed analyses. Figure [5] shows that m max ~ 
L°f for p c = 0.96025(20), with D f = 1.850(5), while 
Fig. [6] shows that m max ~ (p — p c )@ in the limit L — > 
oo, with j3 = 0.172(2) (for a plot with higher resolution 
see the supplementary material (SM)). Both exponents 
are clearly different from the values for OP. Indeed, (3 is 
larger than the value 5/36 = 0.1389 for OP, showing that 
the transition is not more abrupt than in OP, as claimed 
in Q, but less so! 

For d = 3 we also studied systems of up to 2 18 sites, 
with roughly the same number of realizations as for 2d, 
and with similar results (see the SM for details): There 
are also important corrections to scaling, if q is taken too 
large, but they decrease strongly when q is taken as small 
as possible. For q = 0.40 we obtain p c = 0.871(1),/? = 
0.51(1),!/ = 0.86(1), and D f = 2.40(1). These values 
satisfy (like the 2-d exponents) the scaling relation Df = 
d—/3/v, and again they are incompatible with OP (where 
p = 0.4170(3), v = 0.8734(5), D f = 2.5226(1) 0). As in 
2-d, P is clearly larger than in OP, indicating that the 
transition is again less sharp, rather than more abrupt. 



In summary, we have shown that coupling two interde- 
pendent networks does not generically make the perco- 
lation transition more abrupt or discontinuous. Rather, 
the outcome depends on the network topologies. Real 
networks (e.g. transportation, telephone, ...) often are 
locally embedded in space, thus their behavior might re- 
semble more that of regular lattices than that of small 
world networks. The reason why the claim of [H does not 
hold universally is not that the cascade picture breaks 
down for local networks. Rather, cascades are an es- 
sential ingredient in any spreading phenomena on any 
network, and it depends on the topology whether or not 
their effects are enhanced by the coupling between differ- 
ent networks. 

In the present paper we have only studied two statis- 
tically identical networks. It is an open question what 
happens, say, when a diluted 2-d lattice is fully coupled 
to an ER network or a scale-free one. Also, one might 
think of more than 2 interdependent networks [lij ]. In 
view of possible applications, one should also study net- 
works that are semi-locally embedded in 2-d space. The 
latter could also be used to study the cross-over from net- 
works with local connections (as in 2-d lattices) to global 
(e.g. ER) networks. A priori, one might expect that there 
exists a tricritical point between these two extremes, or 
that one of them is unstable against even infinitesimal 
perturbations. 
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In this supplementary material, we present several plots. They illustrate claims for which the data are not shown 
in the main paper, or they are plotted in different ways. All necessary information to explain the plots is given in the 
figure captions. 
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FIG. SI. (Color online) Log-log plot of L _1/ '' / dlog(m max }/dp, where the derivative with respect to p is approximated by a finite 
difference quotient with dp = 1/1024. According to the scaling ansatz Eq. (4), this should be a constant at p — p c . It is indeed 
constant (independent of L) at p = 0.96025 within the expected accuracy (we do not show error bars, since the derivative is 
obtained from correlated measurements, making any error estimation difficult), if we assume v = 1.19. Within the quoted error 
bars on Df and /3, this agrees with the hyperscaling relation Df = 2 — P/v. Notice that the OP value v = 4/3 is definitly 
excluded. 
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FIG. S2. (Color online) Log-log plot of (m max )/[(p — p c )^N] versus (p — p c )L < with p c and /? = 0.172 as given in the 
paper and v — 1.19 as obtained from the hyperscaling relation, with Df as given in the paper. If these values and the 
scaling ansatz Eq. (4) were exact, the curves would perfectly collapse onto a single curve that is horizontal in the region 
O(l) < {p — p c )L}/ v < OiL 1 !"}. Deviations mainly result from finite size corrections. The three straight lines indicate the 
slopes we should observe in the power law region, if ft were different from the nominal value ft = 0.172. The full (red) line 
corresponds to ft = 5/36 in OP, and is solidly excluded. The two dashed lines correspond to ±1<t. Notice the very small range 
of the plotted data in the y direction, leading to enormous blow-up of errors and to a much higher significance than in Fig. 6 
of the main text. 
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FIG. S3. (Color online) Critical values p c versus q for coupled 3-d lattices. This plot is to be compared to Fig. 4 of the main 
paper, where results are shown for two dimensions. This time a lower bound on q m j n is given by the value of p c = 0.248812 . . . 
of 3-d bond percolation, while p c for q = 1 equals the critical value 0.31160 ... of 3-d site percolation. As in two dimensions, we 
suggest that the transition is in the same universality class in the entire range q m in < q < 1, while it is in the OP universality 
class for q = 1. 
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FIG. S4. (Color online) Data collapse for (m max )/L D f versus (p — p c )L 1// " for 3-d lattices. As in Fig. 3 of the main paper, each 
set of curves corresponds to one value of q. In contrast to the 2-d case, we needed however slightly different critical exponents 
in order to obtain good collapses for both values of q. More precisely, we used for q = 0.4 the exponent values /3 = 0.51 and 
v = 0.86 quoted in the main paper, while we used /3 = 0.446 and v = 0.77 for q — 0.6. The values of p c used in this plot are 
0.871 (for q = 0.4) and 0.5464 (for q = 0.6). Notice that both sets of curves in the present figure show visible deviations from 
a perfect collapse, but - in contrast to what might be suggested by the figure - we claim that the effective exponents obtained 
for q = 0.4 are much closer to the true ones. This conclusion is mainly based on the following figures (Figs. S6 to S9) which 
would show substantially more scaling violations for q = 0.6 than they do for q = 0.4. We interpret this as a manifestation of 
the cross-over from OP, that should be much more important at q — 0.6 than at q — 0.4. 
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FIG. S5. (Color online) Log-log plot of (m max )/L D ^ versus L for 3-d interdependent percolation at q = 0.40 and fixed values 
of p (analogous plot to Fig. 5 of the main paper) . The fractal dimension is chosen such that the curve for p = p c is horizontal, 
with p c = 0.871. 
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FIG. S6. (Color online) Analogous plot to Fig. SI, but for d = 3 and q = 0.4. We see strong finite size corrections for small 
values of L. For q = 0.6, scaling violations would extend to the largest values of L (see Fig. S7). 
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FIG. S7. (Color online) Same as Fig. S6, but for q — 0.6. The value for v is that used in Fig. S4 to obtain the best data 
collapse. We see now that this good data collapse was spurious. Either v or p c as obtained in Fig. S4 are wrong, but most 
dramatic is the unphysical curve crossings for p < p c . While the consistency of v and p c could be improved at the cost of 
deteriorating somewhat the data collapse in Fig. S4, the curve crossings for p < p c cannot be avoided. 
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FIG. S8. (Color online) Analogous plot to Fig. S2, but for d = 3 and q = 0.4. The values of p c and of the critical exponents 
are those given in the caption to Fig. S4 and in the main text. 




FIG. S9. (Color online) Analogous plot to Fig. S8, but for q — 0.6. The values of p c and of the critical exponents are those 
given in the caption to Fig. S4 and used also in Fig. S7. The straight line represents the slope expected for OP (since we do 
not assign error bars to the value of /3 used in this plot, we do not show the two other straight lines that were shown in Figs. 
S2 and S8). Again we see that the data collapse in Fig. S4 for q = 0.6 is spurious (it is restricted to (p — p^L 1 ^" < 1), and that 
corrections to scaling are huge at q — 0.6. Similar but less dramatic scaling violations were also seen in d — 2. if q = 0.75 was 
used instead of q = 0.6 (data not shown) . This supports our conclusion that scaling violations are due to cross-over from OP. 



